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Abstract 

The pattern of isolation-by-distance arises when population differentiation increases with increasing 
geographic distances. This pattern is usually caused by local spatial dispersal which explains why dif- 
ferences of allele frequencies between populations accumulate with distance. However, the pattern of 
isolation-by-distance can mask complex variations of demographic parameters. Spatial variations of de- 
mographic parameters such as migration rate or population density generate non-stationary patterns of 
isolation-by-distance where the rate at which genetic differentiation accumulates varies across space. Bar- 
riers to gene flow are particularly well studied examples that generate non-stationary patterns of isolation- 
by-distance. Using the concept of genetic friction, we develop a statistical method that characterizes 
non-stationary patterns of isolation-by-distance. Genetic friction at a sampled site corresponds to the 
local genetic differentiation between the sampled population and Active populations living in the neigh- 
borhood of the sampling site. To avoid defining populations in advance, the method can also be ap- 
plied at the scale of individuals. The inference of genetic friction relies on a matrix of pairwise sim- 
ilarity between populations or individuals. Typical matrices of similarity include correlation matrices 
or i<5T-based matrices. The proposed framework is appropriate for dealing with massive data because 
it relies on a pairwise similarity matrix, which can be obtained with computationally efficient meth- 
ods. A simulation study shows that maps of genetic friction can detect barriers to gene flow but also 
other patterns such as continuous variations of gene flow across habitat. The potential of the method 
is illustrated with 2 data sets: genome-wide SNP data for the human Swedish populations, and AFLP 
markers for alpine plant species. The software FRICTION implementing the method is available at 
|http : / /membres-timc . imaq . fr /Michael . Blum/Friction . html[ 
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INTRODUCTION 



Characterizing the pattern of genetic differentiation within a species is a recurring task in population genetics analyses. 



WRIGHT ( 1943 1 introduced the model of isolation by distance (IBD) which assumes that differences of allele fre- 
quencies between populations accumulate under the assumption of local spatial dispersal. Because of local dispersal, 
IBD models predict the so-called pattern of IBD where population differentiation increases with increasing geographic 
distances ( |SLATKIN|1993||R0USSET|1997| ). This pattern is observed in many model and non-model organisms as well 
as in humans suggesting that local dispersal is a leading evolutionary force (Shar bel et al. |2000[ |RAMACHANDRAN| 
et al. |2005 ||Hardy et aZ.|2006l|HELLBERG|2009| >. 

However, the pattern of IBD can mask complex variations of demographic parameters resulting in differential in- 
creases of genetic differentiation in different regions of the habitat. Variations of demographic parameters can arise 
when population densities or migration rates vary across space ( |Slatkin|1985[ ). With the advent of landscape ecology 
( |PlCKETT and Cadenasso|1995] > and landscape genetics ( |Manel et al.\2(X)2>) , the spatial variation of demographic 
parameters is an important topic because spatial heterogeneity (or landscape characteristics) is now recognized to be a 
key factor to explain population differentiation and gene flow ( |McRAE and Be ier 2007 ). Examples of spatial hetero- 
geneity influencing population differentiation include fragmented landscapes in urban and agricultural area where there 



are 'corridors' for gene flow ( Arnaud 2003| |Munshi-SOUTH|2012| i. These habitats can also be characterized by spa 



tially varying local subpopulation size, another factor causing spatial variation of the rate of differentiation between 
populations ( SERROUYA et q/.]|2012 l. Barriers to gene flow, which can be caused by anthropogenic or geographic 



factors, are also emblematic examples of spatial heterogeneity influencing population structure (e.g. CASTELLA et al. 
[20001 [Epps et aZ.|2005||RiLEY et aZ.|2006||GAUFFRE et aZ.|2008||ZArEWSKi e t aZ.|2009| l. Because the identification 
of the barriers to gene flow has attracted considerable attention ( Storfer et al. 2010| l, there is a large variety of statis- 



tical methods dedicated to their detection pSARBUJANl et q/.|19"89l|BQCQUET-APPEL and Bacro|1994[|Dupanloup| 
\et q/.|2002||MANNl et aZ.|2004t|CERCUElL et a/.|2007HCRlDA and Manel|2007||Manel et a/.|2007||SAFNER et al.\ 
201 Here, we propose a more general method that characterizes non-stationary patterns of isolation-by-distance. 



A non-stationary pattern of isolation-by-distance occurs when the spatial rate of differentiation between individuals 
or between populations depends on space. Non-stationary patterns of IBD arise when there is a barrier to gene flow 
because there is a larger rate of genetic differentiation in the region of the barrier but also on different situations such 
as continuous variations of gene flow across habitat. 

To characterize non-stationary patterns of IBD, our approach provides a measure of local differentiation at each 
location where genetic data is available. By analogy with the concept of friction map in ecology ( |Ray|2005) , which 
measures the cost of movement through the landscape, the measure of local differentiation can be thought of as a 
measure of genetic friction with larger friction values indicating larger population genetic differentiation per unit of 
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spatial distance. The principle of the method is to estimate for each sampled location Zj, i = 1, . . . ,n, a pairwise 
measure of population differentiation or of dissimilarity between the population located at the sampled location and 
Active populations located in the neighborhood of z. t (see Figure[T]i. The method is not restricted to pairwise population 
measurements and can also accommodate individual pairwise measures. Working at the scale of individuals is a 
desirable feature since using individuals as the operational unit avoids potential bias in identifying populations in 
advance and offers the opportunity to conduct studies at a finer scale ( |Manel et q/.| |2003 , 2007 ). Using a detailed 
simulation study, we demonstrate that the method can correctly infer local variation of genetic differentiation and we 
present applications of the method to genome- wide SNP data for the human Swedish population (HUMPH REYS et ai\ 
|2011| l and AFLP makers from alpine plants ( |GUGERLI et aZ.|2 008 ). 



METHODS 

We assume that the data consist of the correlation matrix of allele frequencies over populations. To assess local genetic 
differentiation around a given sampled site, we estimate the correlation of allele frequencies between the sampled 
population and the neighboring populations that have not been sampled. Neighboring populations are located at a fixed 
and short distance from the sampled populations, and we measure the expected correlation (averaged over neighbors) of 
allele frequencies between the sampled population and the neighboring Active populations. Since we aim at providing 
friction values that should be larger in regions of abrupt genetic changes, we consider one minus these local correlations 
as friction measures. 

We estimate the local correlation using kriging, which refers to geostatistical techniques to interpolate the value of 
a variable at unsampled locations from observations of its value at nearby and sampled locations ( Cressie|1993| |. We 
do not use kriging in a standard manner, which would consist of interpolating the allele frequencies at the neighboring 
sites where the Active populations are located. We instead provide an estimation of the correlation of a variable (allele 
frequency here) between the sampled sites and the neighboring sites. The kriging procedure is described below. 

The kriging approach: A viewpoint of kriging is to model the joint values of the variable at sampled and unsam- 
pled sites as a multivariate Gaussian variable ( BlSHOP|2006 i 



{X,Y) ~* M{m,V), (1) 



where 

/ 

* = 



^xx ^xy 



where X (resp. Y) is the vector of values at the sampled (resp. unsampled) sites, to is a constant mean over the habitat 
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(ordinary kriging), ^ xx (resp. ^ yy ) is the covariance matrix of X (resp. Y) and ^ xy contain the covariances between 
the elements of X and Y . The interpolation of the variable at the unsampled sites is obtained using the conditional 
distribution of Y given X, which can be written in the following regression form 

Y — m = t(X - m) + e, (2) 



where t = ^ xy ^ x ^ and e is a residual independent of X (Le and ZlDEK 1992 1. Using the regression equation ^ that 



relates the variable values at sampled and unsampled sites, we can show that the the predictive value of the covariance 
matrix between X and Y is given by 

^->xy — 7~^i xx . (3) 

where H xx is the covariance matrix for the sampled sites and we consider the empirical covariance matrix S in our 
computations. Providing the correlation instead of the covariance between sampled and unsampled sites requires 
the standardization of the covariance equation ([3]) and the renormalization formula is provided in Appendix A. The 
predictive distribution of the matrix Yj xy is finally obtained by integrating over the posterior distribution of r, which 
is a function of the matrix ^. The parametric model for ^ and the sampling algorithm that provides its posterior 
distribution is given below. 

A model for the correlogram: 

We consider the standard model of stationary kriging which assumes that the correlation between two points only 
depend on the distance between these two points. Using these assumptions, we should model how the correlation 
decreases with increasing distance. We assume that this function C, called the correlogram, decays exponentially 

C(d) = (a + (1 - a)e- d / r + Al d=0 )/(1 + A), (4) 

where d is the distance between two points, 1 denotes the indicator function, a determines the sill, which measures the 
limiting value of the correlation function, r is the range parameter and A is the regularization parameter. The parameter 
A is introduced to account for the uncertainty of measure at the sampled sites (BISHOP 2006). More importantly, it 



ensures that the matrix $ xl is invertible, which is required for computing the regression matrix r. The range parameter 
r is inversely related to the rate at which correlation decays with distance. We sample (a, A, r) from the posterior 
distribution using a Gibbs sampling algorithm with details provided in Appendix B ( |HANDCOCK and Stein|1 993 ). 
It might appear paradoxical to consider the model of stationary kriging to provide a spatial estimation of the different 
rates at which correlation decays with distance. However, the stationary model of equation |4]) only provides values 
for the regression matrix r and it is the data contained in Y, xx that convey information about non-stationarity through 
equation ([3]). 
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SIMULATION STUDY 



In the simulation study, we consider two different models for generating non-stationary patterns of isolation-by- 
distance. First, we consider non-homogeneous stepping stone models in one and two dimensions. We simulate with 
ms ( Hudson]20 02) spatially-dependent effective migration rate AN m where iV is the population size of each deme 
and m is the migration rate per generation between two neighboring demes. The second model is analytic and has 
been developed for performing non-stationary kriging when the correlogram function (equation^) is assumed to vary 
across space ( |PAClOREK and SCHERVISH|200 6). The range parameter r of equation Q, which measures the rate at 
which correlation decays with distance, is assumed to be a function of space. For the second model, zone of abrupt 
changes such as genetic barriers correspond to regions with a smaller range parameter because correlation decays more 
rapidly with distance in these regions. 

Barrier in one-dimensional habitat: We provide an example of a one dimensional barrier to show that the param- 
eters of the correlogram function (equation affect the estimated values of genetic friction. We simulate a stepping 
stone model with 100 populations of effective sizes No = 1000 diploid individuals. We sample 20 equidistant popula- 
tions and we consider 20 chromosomes in each of them. Migrations are constant between neighboring populations and 
ANqtu = 20. The barrier is located between populations 50 and 51 and arose 8 units of time ago (4^0 m = 0) where 
time is counted in units of 4A^o generations. As input data (E xx with n = 10), we consider the pairwise correlation of 
allele frequencies for the 20 sampled populations. 

Figure[2]displays the estimated friction values along the one-dimensional habitat when sampling 20 equally-spaced 
populations. For each sampled deme, the friction corresponds to one minus the expected correlation between the sam- 
pled deme and its two neighbors. Whatever the values of the correlogram parameters, (a, A, r), the friction values are 
larger in the middle of the habitat, which is consistent with the presence of a barrier to gene flow. Nonetheless the 
detailed trajectory of the friction curve depends on the choice of the parameter value. To account for the uncertainty 
associated with the parameters of the correlogram function, we integrate the friction values over the posterior distribu- 
tion of the correlogram model (thick line in Figure [2]). When considering a random and uniform sampling of the 20 
populations instead of a regular sampling, we also find that the friction values are larger around the barrier to gene flow 
(Figure SI). 

Barriers in two-dimensional habitat: We consider two examples of a 2-dimensional habitat with genetic barriers. 
In the first example of a 2-dimensional habitat, the data are simulated using a stepping-stone model with a 10 x 10 
grid. In each population, there are N a = 1000 diploids individuals per population and we sample all populations 
considering 20 chromosomes in each of them. The migration rate between neighboring populations is of 4N m = 20 
where migration only occurs along horizontal and vertical lines but not along diagonals. We then assume that two 
barriers arose 7\ = 5 and T 2 = 3 units of time ago where time is counted in units of generations (see Figure|3]l. In 
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the second 2-dimensional example, we specify explicitly the local decay of correlation using the non-stationary model 
of PACIOREK and SCHERVISH (2006). We assume that there are three genetic barriers, which correspond to three 
different regions where the range parameter (r in equation^) is smaller (Figure|4]i. In both examples, the input matrix 
of pairwise values (E xx ) is the correlation matrix for the sampled sites, although the correlation is estimated using 
the simulated allele frequencies for the stepping stone example and it is obtained analytically using the convolution 
formula of |PAClOREK and SCHERVISH| ( |2006| > for the second example. 

Figures [3] and [4] show that the friction values are larger around the genetic barriers as expected. For both examples, 
the relative importance of the barriers is retrieved. The strongest barrier has the largest friction value and the weakest 
barrier has the smallest friction values. For computing friction values, we also consider Fst pairwise values instead of 
correlation values in the stepping-stone model. The matrix of input data T, xx contains the pairwise (1 — Fst) values 
that decay with increasing geographical distance as assumed by the correlogram equation The friction values now 
correspond to the expected Fst between the sampled populations and their Active neighbors. We find that the friction 
map obtained with the Fst measures is similar to the friction map obtained with the pairwise correlation (Figure S2). 

A gradient of gene flow: We also consider a different pattern of non-stationary IBD consisting of a 2 dimensional 
stepping-stone model with a spatial gradient of gene flow. We assume that gene flow is maximum at the lower left 
corner of the habitat and decreases proportionally with distance from the lower left corner of the habitat (Figure[5]l. As 
expected, we find that genetic friction increases when moving away from the lower left corner of the habitat (Figure[5]l. 

This example is particularly interesting instance of non stationarity, which can not be described with barriers to 



gene flow. For the previous examples, the software barrier, which detects zones of abrupt genetic change ( Manni 



et al. 



2004), was able to retrieve the barriers in both the one and two-dimensional habitat (Figure S3). However, for the 



gradient of gene flow, barrier incorrectly finds a barrier in the upper right corner of the habitat, which is nonetheless 
consistent with the fact that gene flow is minimal here (Figure [5]). Additionally, multidimensional scaling — another 
commonly used method to represent differentiation between populations — provides a meaningful representation for 
the examples of barriers because the populations that live on different sides of the barriers are clearly separated in the 
two-dimensional representation of MDS (Figure S4). However, the interpretation of the pattern obtained with MDS is 
much more difficult for the example of a gradient of gene flow (Figure [3J. The observed pattern found with MDS is 
consistent with the gradient of gene flow because populations living in regions of high gene flow (dark points in Figure 
[5]) are located more closely on the MDS plot than populations living in regions of low gene flow (clear points in Figure 
[5]). Although consistent with a gradient of gene flow, the MDS plot is not as easily interpretable as the map of genetic 
friction for this example. 
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APPLICATIONS 



Non-stationary patterns of IBD among the Swedish population: We first illustrate the kriging methodology using a 
human SNP data set with a particularly dense geographic sampling. The data consist of genome-wide SNPs for 5174 
Swedish individuals that cover all of the 21 Swedish counties (HUMPHR EYS et a/.|201 1 1. To assign each individual to 
a county, HUMPHRE YS et q/.| ( |201 l| l used available geographic information with the following order of priority: city or 
village of birth, county of birth, municipality or city of residence and county of residence if it is the only information 
available. They found strong differences between the far northern counties and the remaining counties, and additionally 
demonstrated that the northern counties are more clearly genetically differentiated from each other than the southern 
counties are from each other. 

Since our framework is an extension of isolation-by-distance, we first check that population differentiation increases 
with increasing geographical distance. We confirm the prevalence of isolation by distance in Sweden (P < 10~ 7 for 
a Mantel test, see also Figure S5). Then, we choose to quantify genetic friction using the F$t between a population 
living exactly in the barycentric center of the county and a putative neighboring population living 30 km away (see 
Figure S6). We consider the pairwise (1 — F$t) values between the counties as input matrix of pairwise similarities 
(T, xx ). We find that the northernmost counties (Nordbotten, Vasterbotten and Jamtlands) have the strongest friction 
values whereas the smallest values were found in the regions around the Stockholm area (Ostergotlands, Stockholms, 
Stidermanlands, Vastmanlands and Jonkopings are the five counties with the lowest friction values). The fourth largest 
value is in the Dalarna county (in north middle Sweden), which borders southern Norway, and counts more individuals 
with remote Finnish or Norwegian ancestry than other counties ( HUMPHREYS et aT||2011 1. As expected, the four 
counties with the largest friction values are also the most differentiated from other counties even when controlling for 
geographic distance (Figure S5). 



In summary, we confirm the results of HUMPHREYS et al. (201 1 1 who found that there is more genetic differenti- 
ation within northern Sweden than within southern Sweden. The four counties with the largest friction values are also 
the counties with the lowest population densities (Figure S7) suggesting that low density habitat triggered population 
differentiation in Northern Sweden. 

Non-stationary patterns of IBD for 20 alpine plants: We consider a set of 20 alpine plant species that have been 
sampled across the Alps ( |GUGERLI et a/.|[2008{ | ALVAREZ et a/.||2009) |JAY et a/.||2012[ l. The sampling is particu- 
larly dense with three individuals per species collected and genotyped for each cell of approximately 500 km 2 . The 
individual genotype consists of amplified fragment length polymorphisms (AFLPs). Because of the individual-based 
sampling, we consider the matrix of correlation between the individual genotypes as input data matrix (E xx ). The 
genetic friction that we estimate corresponds to the expected correlation between sampled individuals and neighboring 
individuals located at 10 km. 
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The map of genetic friction for Rhododendron ferrugineum, one of the 20 sampled plant species, exhibits large 
friction in the southern part of a corridor that cuts the main chain of the Alps and roughly stretches from Innsbruck 
to Lake Garda (Figure |7J. This pattern is representative of what is generally found for the other species because this 
region corresponds to the main zone of genetic friction for 19 out of the 20 alpine species we consider (see Figure S8 
for all friction maps). When using a log scale for the friction values, we find for most of the species that the whole 
corridor, and not its Southern part only, has more elevated friction values than the rest of the Alps (Figure S9). The 
southern part of the corridor corresponds to the Adige valley, which is perpendicular to the main E-W axis of the Alps 
and which further leads to the Reschen and the Brenner pass, the later being the lower pass across the main chain of 
the Alps ( |TAYLOR|1940| l. For plant species, the Adige valley — also referred as the Brenner line — is already known 
to be an important barrier for plants both at the intraspecific and interspecific level ( Thiel-Egenter et a/.||2011| l. 
The explanation for the presence of a barrier in the Brenner line involves Pleistocene glaciations: the populations of 
plants were initially fragmented into glacial refuges, then expanded via postglacial colonization routes, and a secondary 
contact zone finally arose around the Brenner line where formerly allopatric populations admixed ( |PAWLOWSKl|1970| 
ISCHONSWETTER et a/.|2005[[T~HIEL-EGENTER et a/.|2011) . 

To check that the observed pattern is not an artifact of the sampling scheme, we simulate a pairwise matrix of 
correlation values (E xx ) using the stationary model of equation Q which assumes that correlation decays at the same 
rate in the entire habitat. We consider values for the parameters of the correlogram that provide a good fit to the 
observed decay of correlation in Rhododendron ferrugineum. The inferred map of genetic friction indicates that the 
kriging method correctly infers homogeneous friction, which confirms that the large friction found around the Adige 
river is not an artifact of the kriging method (Figure [7}. 

DISCUSSION 

In this article we present a new Bayesian method to characterize non-stationary patterns of isolation-by-distance. From 
global measures of pairwise similarity or dissimilarity, the method infers local measures of similarity or dissimilar- 
ity. Whatever is the exact measure of genetic (dis)similarity, we use the generic expression of genetic friction when 
referring to the estimated local growth of genetic dissimilarity or differentiation. If considering for instance the F$t 
pairwise matrix of genetic differentiation between populations, the inferred genetic friction correspond to the F$t 
between the sampled populations and Active neighboring populations located at a given distance. The method is not 
restricted to Fst measures and can handle any type of measures of differentiation and is also valid at the individual 
scale. We consider for instance the correlation between the allelic types of individuals, but other measures would be 
valid such as identity by descent between individuals (BROWNI NG and Browning|2011| i as well as coancestry mea- 
sures ( jLAWSON et a/.|2012| l. Because the two latter measures are based on haplotypes instead of genotypes, they can 
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provide information at a finer scale ( |GATTEPAILLE and JAKOBSSON|2012||LAWSON et al.\20\2) . 

Genetic differentiation and gene flow: It is of course tempting to convert maps of genetic friction into maps 
of gene flow or of dispersal distance. Assuming that differentiation occurs according to a stepping stone model, 
such parameter estimates could be obtained using theoretical relationships between local F$t an d dispersal distance 
(ROUSSET 1997 ). However, it is well known that relating F$t or other measures of genetic differentiation to gene 
flow relies on many assumptions that may be unrealistic ( |MARKO and HART|2011| l. Although the estimation of gene 



flow with F s t-based methods can be robust in some situations, such as temporal variation of gene flow (LEBLOIS 



\et a/.|[2004) l, there are other processes such as range expansion, local extinction and recolonization that can modify 
drastically the pattern of genetic differentiation ( |WADE and McCauley|1988[|ExC0FFIER et qZ.|2009"l|ARENAS et al.\ 
|2012|[SLATKIN and Excoff ier|20T2| i. More generally, providing a map of genetic friction is informative about the 
pattern of genetic differentiation but does not distinguish among the evolutionary processes that generated this pattern 
(for similar concerns about PCA, see | M C V I-: A \ 1 2( )( )9 ) . To provide a concrete example, we find that the same pattern, a 
zone of elevated genetic friction, can correspond to a barrier to gene flow (Figures |2]|4| or to a secondary contact zone 
in the case of the alpine species (Figure|7]i. 

Relevance to landscape genetics: The two key steps of landscape genetics are the detection of genetic discontinu- 



ities and the correlation of these discontinuities with landscape and environmental features such as barriers ( MANEL 



et a/.|2003] l. Detection of genetic discontinuities is clearly provided by the Kriging method; for instance in the case 



of the alpine species, we find genetic discontinuities, i.e. larger genetic friction, around an important Alpine valley 
(Figure [7J. The second key step where these discontinuities are correlated with landscape or environmental variable 
can also be obtained as a post-processing step by correlating estimated genetic friction with landscape variables. For 
instance, in the case of the human SNP Swedish data, we find that genetic friction is correlated with population density. 
There are alternative and integrative approaches that account for both genetic data and landscape variables within the 
same statistical framework. Accounting for both sources of data can be performed either by a joint assessment of the 



pattern of population structure or differentiation and its correlation with landscape or environmental variable (FOLL 



and GAGGIOTTI 2006 Jay et a/.|2011| ), or by correlating genetic distances with distances based on landscape features 



( CUSHMAN et a/.|20 06 ; McRae 2006 ). These integrative approaches are hypothesis-driven in the sense that each set 
of landscape features affecting population structure corresponds to one hypothesis that can be tested or compared to 
other ones. The proposed kriging approach is instead a technique of exploratory data analysis. It might be especially 
appropriate for large-scale conservation studies not focused on the underlying evolutionary processes but that should 
deal with reserve design and with the management of fragmented populations (SCHW ARTZ et al.\2Q01\ . To compare 
or test the support of different evolutionary processes provided by the pattern of non-stationary IBD, we should rather 
resort to inference based on explicit simulations of evolutionary processes using for instance approximate Bayesian 
computation ( CSILLERY et a/. |2010 1. Within this simulation framework, measures of genetic friction can be included 
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as statistical summaries of the data. 

Caveats: As with other descriptive techniques of exploratory data analysis, there is no hypothesis-testing or good- 
ness of fit procedures attached to it. This is a concern of primary importance because it can lead to overinterpretations 
of friction maps. When analyzing the alpine plant species data, we proposed one solution. We simulated spatially 
homogeneous patterns of genetic differentiation using the same sampling scheme as in the data and we checked that 
the obtained friction pattern is indeed more homogeneous than the one obtained with the data (Figure [7). Using the 
same sampling data as in the data is crucial since the sampling scheme affects the ascertainment of population structure 
for various statistical methods ( |McVean"||2009[ |SCHWARTZ and McKelvey|[2009| >. An additional concern about 
the kriging method concerns the choice of the distance between the sampled populations or individuals and the Ac- 
tive neighboring populations or individuals. We considered various choices of distances for both the Swedish data set 
(10-150 km) and one of the alpine species (5-100km) and found that these large ranges of distances provided similar 
patterns of genetic friction although the patterns generated with the larger distances were smoother (Figures S10 and 
Sll). The last caveat to bear in mind concerns the choice of the pairwise dissimilarity or differentiation matrix. For 



instance, in a simulation study with 5 populations, LAWSON and FALUSH (2012 i showed that the five populations were 
clearly distinguishable with some but not all pairwise dissimilarity matrices between individuals. Although a poten- 
tial caveat, being able to choose the measure of dissimilarity also adds to the flexibility of the method and different 
measures can convey information about processes that occurred at different time periods. 

Perspectives: When carefully addressing the aforementioned caveats, measures of genetic friction can provide in- 
terpretable patterns for describing non-stationary patterns of IBD. The method relies on a matrix of pairwise (dis)similarity 
between individuals or populations located on georeferenced sampling sites. An approach based on dissimilarity ma- 
trices is an appropriate methodology for the new genomic era where we have to deal with massive data. Computing 
pairwise dissimilarity matrix can be computationally efficient (e.g. BROWNING and Browning |2011 i and can even 
be parallelized to compute different parts of the matrix ( LAWSON and Falush|[20T2| . Having a statistical method 
able to scale with the dimension of the genetic data should make it a valuable tool for investigating patterns of genetic 
differentiation in a wide range of studies. 
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APPENDIX A: COMPUTING THE CORRELATION BETWEEN SAMPLED AND 

UNSAMPLED SITES 



To renormalize equation we have to compute the diagonal elements of the covariance matrix T, yy for the unsampled 
sites. Using the fact that the residuals of the regression equation (|2j are of variance ^ yy — ^ xy^xx^ xy (BISHOP 2006 1, 
we can show that the covariance matrix T, yy is given by 



^yy '"xx' <^yy *%y*xx*xy 



APPENDIX B: GIBBS SAMPLER 

For the parameters of the correlogram model, we choose the following prior distributions 

a ~ hi (min — <5, min + 5) (5) 

A ~ log(W(0, 1)) (6) 

r ~ U (min Z3y , max Dij ) , (7) 

where S is the empirical covariance matrix and D denotes the pairwise geographical distance between sampled sites. 
Here, we adopt an empirical Bayes approach, since we partly use the data when defining the prior distributions. 
Using the Bayes formula we have 

log(P(a, A, r\S)) = log(p(S\a, A, r)) + log(p(a, A, r)) + C x , (8) 



where C\ is a constant and the likelihood is given by the Wishart distribution (SCHMI DT and O'HAGAN|2003| l 



]og(p(S\a, A, r)) - -^log(|*«*l) - \^e(^S) + C 2 , (9) 

where C2 is an other constant. We then simulate a sample of the joint posterior probability using the Gibbs sampling 
algorithm. To obtain replicates from each conditional distribution, we consider a grid for each parameter. We evaluate 
the conditional densities for each point of the grid using equations <JSJ and Q and then we use a sampling algorithm 
for simulating discrete random variables with known probability masses. 
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Figure 1 : A 2-dimensional habitat with putative sampled locations and neighboring locations in grey. Local differenti- 
ation at each sampled location corresponds to the average (over neighbors) pairwise measure of differentiation between 
the population or individual at the sampled location and its neighbors. 
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Figure 2: Estimation of the friction in a one-dimensional habitat with a genetic barrier in the middle of the habi- 
tat. Genetic friction corresponds to one minus the correlation of allele frequencies between sampled and unsampled 
neighboring demes. The genetic barrier is located between the deme 50 and 51. 
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Figure 4: True and estimated friction maps in a two dimensional habitat of dimension 15 x 10 containing three barriers. 
The friction at sampled points corresponds to one minus the correlation between sampled points (large crosses in 
panel A) and Active neighbors (small crosses in panel B) located at a distance of 0.1. The correlation landscape 
is specified with the non-stationary model of |PAC!OREK and SCHERVISH| (2006) considering smaller range values 
(parameter r of equation |4]i for the zones of abrupt genetic change. The true friction of panel A is obtained using the 
convolution formula of PACIOREK and SCHERVISh] ([2006 ) at distance 0.1 which provides pairwise correlation based 
on the landscape of the range parameter of equation d4L 
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A. Friction map and barrier 



B. Multidimensional scaling 
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Figure 5: Investigating the pattern of population differentiation for a gradient of gene flow in a two dimensional habitat. 
Gene flow is maximal at the lower left corner of the habitat and decreases proportionally to the distance from the lower 
left corner. Panel A) Estimation of friction values using the correlation with neighbors located at a distance of 0.1. The 
first genetic barrier that is found with the software barrier is shown with a thick black line. Panel B) Multidimensional 
scaling plot with a grayscale color scheme to represent the distance of each population to the lower left corner. Dark 
points correspond to populations that are close to the lower left corner whereas the lightest points are the farthest away 
from the lower left corner. 



22 




Figure 6: Map of genetic friction computed for the 20 Swedish counties. Genetic friction corresponds to the F st 
between sampled county and Active neighboring populations, not shown, located at 30 km. 
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Data simulated with a stationary process 
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Figure 7: Map of genetic friction computed across the Alps for Rhododendron ferrugineum. Genetic frictions corre- 
sponds to the correlation between sampled individuals and Active neighboring individuals, not shown, located at 10 km 
around the sampled individuals. 
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